1. Root Finding#

강좌: 수치해석 프로젝트

이론#

비선형 방정식의 해를 찾는 과정으로 반복적으로 해석한다.

대표적인 방법은 Bisection Method, Newton Rapson Method가 있다.

Bisection Method#

  • 해가 구간 (a, b) 사이에 있을 때 해가 존재하면 \(f(a)\cdot f(b) < 0\) 조건을 활용한다.

  • 이 경우 중간 값 \(c=\frac{a+b}{2}\) 에 대해

    • \(f(a)\cdot f(c) < 0\) 이면 해가 (a, c) 사이에 있다고 범위를 좁힘

    • \(f(c)\cdot f(b) < 0\) 이면 해가 (c, b) 사이에 있다고 범위를 좁힘

Note

\(|f(c)| < tol\) 이면 멈춘다.

bisect-fig

Fig. 1 Bisection method (From Wikipedia)#

Python 구현#

def bisect(a, b, f, tol=1e-12):
    """
    Bisection method
    
    Parameters
    ----------
    a : float
        Lower limit
    b : float
        Upper limit
    f : function
        함수
    tol : float
        Tolerance
    """
    product = f(a)*f(b)
    
    if product > 0:
        # 같은 부호 이므로 이 구간 내에 해가 존재하지 않음
        raise ValueError('Wrong Intervals')
    else:
        c = 0.5*(a+b)
        
        if abs(f(c)) < tol:
            print('Converged at {:.7g}'.format(c))
        elif f(a)*f(c) < 0:
            # 해는 (a,c) 사이에 존재함, 범위 좁히기
            bisect(a, c, f, tol)
        else:
            # 해는 (c,a) 사이에 존재함, 범위 좁히기
            bisect(c, b, f, tol)

예제#

%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np

plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
x = np.arange(-10, 10, 0.1)

def f(x):
    return x**2 + 10*np.sin(x)
plt.plot(x, f(x))
[<matplotlib.lines.Line2D at 0x7f49dc146210>]
../_images/5405813438f39ea3dfae2e33ce0a05626e30dac8e12faab0d516c95e036a1879.png
# Left solution
bisect(-10, -1, f)
Converged at -2.479482
# Right solution
bisect(-1, 1, f)
Converged at 0
# Wrong Interval 1
bisect(-10, -5, f)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [7], line 2
      1 # Wrong Interval 1
----> 2 bisect(-10, -5, f)

Cell In [1], line 20, in bisect(a, b, f, tol)
     16 product = f(a)*f(b)
     18 if product > 0:
     19     # 같은 부호 이므로 이 구간 내에 해가 존재하지 않음
---> 20     raise ValueError('Wrong Intervals')
     21 else:
     22     c = 0.5*(a+b)

ValueError: Wrong Intervals
# Wrong Interval 2
bisect(-5, 5, f)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [8], line 2
      1 # Wrong Interval 2
----> 2 bisect(-5, 5, f)

Cell In [1], line 20, in bisect(a, b, f, tol)
     16 product = f(a)*f(b)
     18 if product > 0:
     19     # 같은 부호 이므로 이 구간 내에 해가 존재하지 않음
---> 20     raise ValueError('Wrong Intervals')
     21 else:
     22     c = 0.5*(a+b)

ValueError: Wrong Intervals

Newton Rapson Method#

  • 함수 \(f(x)\) 가 미분 가능하고 연속함수인 경우에 대해서 다음과 같은 근사식을 만족한다.

\[ 0 = f(p_0) + f'(p_0)(p_1 - p_0). \]

즉,

\[ p_1 = p_0 - \frac{f(p_0)}{f'(p_0)} \]

이 과정을 반복하면 해를 구할 수 있다.

\[ p_n = p_{n-1} - \frac{f(p_{n-1})}{f'(p_{n-1})} \]

Note

\(|f(p_n)| < tol\) 이면 멈춘다.

  • 초기 Guess 값에서 시작해야 하며, 이에 따라 해가 달라진다.

  • 방정식의 미분 \(f'(x)\)가 필요하며, 이론적으로 구하기 어려운 경우 수치적으로 계산해야 한다.

\[ f'(x) \approx \frac{f(x+\epsilon) - f(x-\epsilon)}{2\epsilon} \]
newton-fig

Fig. 2 Newthon Method (from Wikipedia)#

Python 구현#

def newton(f, df, x0, tol=1e-9):
    """
    Newton Rapson method
    
    Parameters
    ----------
    f : function
        함수
    df : function
        도함수
    x0 : float
        초기 Guess
    tol : float
        Tolerance
    """   
    # Do It Yourself
    if abs(f(x0)) < tol:
        print('Converged at {:.7g}'.format(x0))
    else:
        # x 절편 (Intercept)를 계산함
        x0 -= f(x0) / df(x0)
        
        # 새로운 점에 대해 Newton Rapson 기법 적용
        newton(f, df, x0, tol)

예제 해석#

도함수는 이론적으로 또는 수치적으로 계산할 수 있다.

def df_exact(x):
    return 2*x + 10*np.cos(x)


def df_apprx(x, tol=1e-10):
    return (f(x+tol) - f(x-tol)) / (2*tol)
# Left solution
newton(f, df_exact, x0=-4)
Converged at -2.479482
# Left solution with
newton(f, df_apprx, x0=-4)
Converged at -2.479482
# Right solution
newton(f, df_exact, x0=3)
Converged at 5.448369e-16

Scipy 활용#

scipy.optimize 모듈은 최소화, Curve fitting 그리고 root finding과 관련된 다양한 알고리듬을 제공한다.

참고

from scipy import optimize

scipy.optimize.root_scalar 또는 scipy.optimize.root 함수를 이용하면 손쉽게 계산할 수 있다.

Note

scipy.optimize.fsolve 함수는 MINPACK 라이브러리의 wrapper로 MATLAB의 fsolve와 용법이 같으나, 이제는 안 쓸 것을 권장한다.

# Root with Hybrid
optimize.root(f, x0=-4)
 message: The solution converged.
 success: True
  status: 1
     fun: [ 1.510e-14]
       x: [-2.479e+00]
    nfev: 8
    fjac: [[-1.000e+00]]
       r: [ 1.285e+01]
     qtf: [ 1.069e-08]
# Root_scalar with brentq method
optimize.root_scalar(f, bracket=[-10, -1])
      converged: True
           flag: 'converged'
 function_calls: 12
     iterations: 11
           root: -2.479481833541344
# Root_scalar with Newton Rapson method
optimize.root_scalar(f, x0=-4, fprime=df_apprx)
      converged: True
           flag: 'converged'
 function_calls: 10
     iterations: 5
           root: -2.479481833541344
# Find minimum from x0=0
optimize.minimize(f, x0=0)
  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: -7.945823375615215
        x: [-1.306e+00]
      nit: 5
      jac: [-1.192e-06]
 hess_inv: [[ 8.589e-02]]
     nfev: 12
     njev: 6
# Find minimum from x0=5
optimize.minimize(f, x0=5)
  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: 8.31558557947746
        x: [ 3.837e+00]
      nit: 5
      jac: [ 2.384e-07]
 hess_inv: [[ 1.188e-01]]
     nfev: 12
     njev: 6